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We propose the general multi-band quasiclassical Eilenberger theory of superconductivity to de¬ 
scribe quasiparticle excitations in inhomogeneous systems. With the use of low-energy projection 
matrix, the M-band quasiclassical Eilenberger equations are systematically obtained from A-band 
Gor’kov equations. Here M is the internal degrees of freedom in the bands crossing the Fermi 
energy and N is the degree of freedom in a model. Our framework naturally includes inter-band 
off-diagonal elements of Green’s functions, which have usually been neglected in previous multi¬ 
band quasiclassical frameworks. The resultant multi-band Eilenberger and Andreev equations are 
similar to the single-band ones, except for multi-band effects. The multi-band effects can exhibit 
the non-locality and the anisotropy in the mapped systems. Our framework can be applied to an 
arbitrary Hamiltonian (e.g. a tight-binding Hamiltonian derived by the first-principle calculation). 
As examples, we use our framework in various kinds of systems, such as noncentrosymmetric super¬ 
conductor CePtsSi, three-orbital model for Sr 2 Ru 04 , heavy fermion CeCoIns/YbCoIns superlattice, 
a topological superconductor with the strong spin-orbit coupling CuzBhSes, and a surface system 
on a topological insulator. 

PACS numbers: 


I. INTRODUCTION 

Multi-band superconductors such as MgB 2 and the 
iron pnictides have been attracted much attention be¬ 
cause of its high critical temperature. Although MgB 2 is 
a phonon-mediated superconductor, it has a large critical 
temperature T c ~ 40K, originating from the multi-band 
effects. Multi-band effects are recognized as one of the 
ways to increase the critical temperature. The discov¬ 
ery of the iron-pnictides had a striking impact on many 
researchers in condensed matter physics. Many kinds of 
phenomena with multi-band effects have been proposed 
and confirmed in iron-based superconductors |ljj3j]. Fur¬ 
thermore, recently found superconductors characterized 
by topological invariants, i.e. topological superconduc¬ 
tors, are also multi-band superconductors, since the in¬ 
ternal degrees of freedom {e.g. spins, orbitals or particle- 
hole spaces) in a multi-band system induce topological 
twists in wave functions0-0|. 

The huge computational cost originating from the mul¬ 
tiple degrees of freedom prevents theorists from under¬ 
standing the physical properties in multi-band supercon¬ 
ductors. For example, in the iron-based superconductor 
LaFeAsO, the five-orbital two-dimensional tight-binding 
model has been used as the effective model j3|]. There 
are also ten-orbital three-dimensional tight-binding mod¬ 
els as effective models to analyze experiments in another 
iron-pnictides. {13 In addition, when dealing with vor¬ 
tices and surfaces in multi-band systems, the computa¬ 
tional cost becomes huger, since the momentum is not a 
good quantum number in inhomogeneous systems. For 
example, in the topological superconductors, it is im¬ 
portant to study the quasiparticle excitations so-called 
the Majorana fermions around surfaces and vortices, in 
terms of the bulk-edge correspondence [|j . The Ginzburg- 
Landau framework, which is usually used to examine 


the distribution of the order parameter in the inhomo¬ 
geneous superconductors, is not suitable for dealing with 
the quasiparticle excitations. Even if we use the mean- 
field framework such as the Bogoliubov-de Gennes frame¬ 
work, the simulations of the nano-size multi-band super¬ 
conductors needs enormous computational costs. 

We point out that the effective models {e.g., derived 
by the first principles calculations) might have too many 
bands to describe the low energy physics of supercon¬ 
ductivity. We should note that the number of the bands 
crossing the Fermi level in normal states is less than four 
even in models for iron-pnictides. The low energy physics 
in superconducting state are characterized by the quasi¬ 
particles on the bands crossing the Fermi level in normal 
states, since a characteristic energy scale of the super¬ 
conducting gap (^meV) is much smaller than that of the 
bands far from the Fermi level (~ eV). In the single-band 
weak-coupling Bardeen-Cooper-Schrieffer (BCS) frame¬ 
work, the theory using information only at the Fermi 
surface, called the quasiclassical Eilenberger theory, has 
many successes [l9| . In multi-band superconductors, 
eliminating the high-energy bands not crossing the Fermi 
level can reduce the number of the bands in a low energy 
effective theory as shown in FigJTJ since the high-energy 
bands can not affect the physical quantities in supercon¬ 
ducting state. 

The quasiclassical Eilenberger theory is successful in 
the BCS model of superconductivity. The theoretical 
framework is based on the fact that the coherence length 
£ is sufficiently greater than the Fermi wavelength 1/fcp- 
Various kinds of analytical and numerical techniques on 
the quasiclassical Eilenberger theory have been developed 
and successfully applied to the studies of a large number 
of conventional and unconventional superconductors [fSj- 
H3. The Eilenberger theory was applied into the two- 
band superconductor MgB 2 . In the conventional models 
for MgB 2 , one neglects the off-diagonal inter-band cle- 
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FIG. 1. (Color online) Schematic figure of the multi-band 
Eilenberger theory. The band dispersions are calculated in 
the five-orbital effective model for LaOFeAsH:, ;f|. The bands 
in the shaded regions are neglected in the multi-band Eilen¬ 
berger theory. 


merits in Green’s function. In this case, two decoupled 
Eilenberger equations can describe the quasiparticle ex¬ 
citations, and the multi-band effects are included only 
through solving the gap equations [2tl . 

There are two kinds of bases to consider the multi¬ 
band systems. First one is the basis which is orthogonal 
in the momentum space (so-called “band”-basis). Other 
is the basis which is orthogonal in the real space, such as 
the d-orbitals in models for iron-based superconductors 
and spins in a model with the spin-orbit coupling term. 
We call these real-space orthogonal basis “orbital” basis 
in this paper. In the quasiclassical Eilenberger theory, 
the quasiclassical Green’s functions depend on the mo¬ 
mentum of the relative motion in momentum space and 
the center-of-mass coordinate in real space. In the pre¬ 
vious multi-band quasiclassical theories [24], j25|, the de¬ 
coupled Eilenberger equations are used by neglecting the 
off-diagonal elements of the Hamiltonian in both band 
basis in momentums space and orbital basis in real space. 
In the case of the two-band model for MgB 2 , this as¬ 
sumption is valid to describe the quasiparticle excita¬ 
tions. There is, however, no Eilenberger theory which 
includes the off-diagonal inter-band elements, except for 
a perturbative approached). The inter-band elements of 
Green’s function are important in the complicated multi¬ 
band systems, such as the iron-pnictides and the topo¬ 
logical superconductors. For example, the iron-based su¬ 
perconductors have many entangled Fe-3d orbitals at the 
Fermi level. The Hamiltonian proposed in the iron-based 
superconductors is diagonalized in momentum space by 
the momentum dependent unitary matrix. The ratio 
of which orbital is dominant originates from this uni¬ 
tary matrix and depends on the momentum. This “or¬ 
bital” character, how the orbitals are entangled, at the 
each Fermi wave number is important to understand 
the physics in iron-based superconductors^- The off- 
diagonal inter-band elements of Green’s functions are in¬ 
duced by those of the unitary matrix. These off-diagonal 
elements become important when a self-energy is induced 


by an inter-band scattering, which is important in a sys¬ 
tem with impurities or vortices. A model of topological 
superconductors usually have the off-diagonal elements 
in spin-space due to a spin-orbit coupling. The spins ro¬ 
tate in momentum space, originating from the spin-orbit 
coupling so that the Hamiltonian is not diagonal with 
the use of spin basis in real space. The “spin” character, 
how the spins are entangled, induces topological super¬ 
conductivity even in system with s-wave on-site pairing 
interaction in two dimension [28l l29j. Therefore, the in¬ 
formation about “orbital” characters in momentum space 
is important factor to describe multi-band effects. 


In this paper, we propose a quasiclassical Eilenberger 
framework in multi-band superconductors with a system¬ 
atic low-energy projection. By eliminating the high en¬ 
ergy bands far from the Fermi level, we derive the multi¬ 
band Andreev equations and quasiclassical Eilenberger 
equations in the projected space constructed from the 
bands crossing the Fermi level. We show that the resul¬ 
tant multi-band Eilenberger equations are similar to the 
single-band ones, except for some corrections to describe 
multi-band effects. The quasiclassical framework uses the 
fact that the coherence length £ is usually longer than the 
Fermi wave length l/k-p in a lot of superconductors. The 
orbital characters on the Fermi surfaces in normal states 
are naturally included in our theory. 


This paper is organized as follows. In Sec. [TTJ we in¬ 
troduce the model of the multi-band superconductors. 
The mean-field multi-band Bogoliubov-de Gennes (BdG) 
Hamiltonian is proposed. We introduce the multi-band 
BdG equations and the gap equations, which is the start¬ 
ing point of the quasiclassical theory. In Sec. Mil we dis¬ 
cuss the decoupled Eilenberger theory used in past stud¬ 
ies. We show that orbital characters can not be included 
in this theory. In Sec. II VI we derive the multi-band quasi¬ 
classical Andreev equations, starting with the multi-band 
BdG Hamiltonian. The Andreev equations describe the 
wave functions in the quasiclassical approach. In Sec. cm 
the multi-band quasiclassical Eilenberger equations are 
derived. The Eilenberger equations are the equations of 
motion of the quasiclassical Green’s function. We dis¬ 
cuss the difference between the previous theory and our 
theory. In Sec. IVII we discuss the physical meanings of 
our multi-band Eilenberger theory. In Sec. IVIII we apply 
the multi-band Eilenberger theory in the various kinds of 
systems as examples. We show that the previous theoret¬ 
ical results are reproduced by our theory and the correc¬ 
tions originating from orbital characters are important in 
multi-band systems. In Sec. IVIII! the summary is given. 
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II. MODEL 


B. Multi-band Gor’kov equations 


A. Multi-band BdG equations 

Let us start with the mean-field BdG Hamiltonian in 
the 27V x 2TV Nambu-Gor’kov space, 

% = i/dnWn^n,^). (l) 

Here, the column vector 'I'(t’) is composed of TV fermionic 
annihilation ip a and creation operators i/d at the posi¬ 
tion r (a = l,--- ,7V), V(r) = ({VgM}> {ipi(r)}) T , 
where {4> a {r)} = ,i/j N {r)) T and {V4( r )} = 

(^{(r), • • • ,4>lj(r)) T . The subscript a in ^ a (r) or ^l(r) 
indicates a quantum index depending on spin or orbital, 
etc. These quantum indices are labeled by the orthogo¬ 
nal basis in real space, which we call orbital basis. The 
Bogoliubov-de Gennes Hamiltonian with a matrix form 
is composed of 


H(r 1 ,r 2 ) = H N (r 1 ,-iV 1 )5(r 1 - r 2 ) + A(r 1 ,r 2 ). (2) 


Throughout the paper, hat a denotes a TV x TV matrix and 
check a denotes a 27V x 27V matrix. The 27V x 27V normal- 
state Hamiltonian matrix H(r i, r 2 ) and superconducting 
order parameter matrix A(r!, r 2 ) are respectively defined 
by 




H N (n,-iVi) 

0 


. 0 


( 3 ) 


A(ri,r 2 ) 


0 A(r 1 ,r 2 )\ 

A t (r 2 ,r 1 ) 0 ) 


( 4 ) 


The Dyson equation in Nambu-Gor’kov space (Gor’kov 
equation) is obtained by 

J dr' (iu n 15(r 1 - r') - H N (r 1: r') 

-A(ri,r*') - G(r', r 2 , 7w n ) = 6(r i - r 2 )l, 

( 8 ) 

with 

iL N (r 1; r 2 ) = (tf N0 (r 1 ) + H m (-iV ri )) <5(n - r 2 ), 

(9) 

where, uj n = (2 n + 1 )ttT is the fermionic matsubara fre¬ 
quency, E(ri,r',7w n ) denotes the self-energy. Here, the 
27V x 27V Green’s function is determined by 

G(r 1 ,r 2 ,T 1 - t 2 ) = — (T T T(ri, ri)'k t (r 2 , r 2 )), (10) 

_ / G(ri,r 2 ,ri - t 2 ) F(ri,r 2 ,n - r 2 ) 
\ F(n,r 2 ,n - r 2 ) G(ri,r 2 ,n - r 2 ) 

(ID 

G Q/3 (ri,r 2 ,n — r 2 ) = i, ri)^(r 2 , t 2 )), (12) 

7A/3(ri,r 2 ,Ti — r 2 ) = —(T T ip a (ri, Ti)ij)p(r 2 ,T 2 j), (13) 

Fa/3(r 1 ,r 2 ,T 1 — t 2 ) = -(T r V’l(ri,ri)^(r 2 ,T 2 )), (14) 

G Q/ 3 (ri,r 2 ,n — r 2 ) = -(T r ^(n, Ti)ipp(r 2 , r 2 )), (15) 

with the imaginary time r. The local density of states is 
given by 

7V(r, 17) = —Im 

7r 

(16) 

III. DECOUPLED EILENBERGER EQUATIONS 


lim Tr G(r,r,iui n ^>E + ir)) 


The order parameter matrix is given by (so-called the 
gap equations) [loj 

A aP {r u r 2 ) = Y, e lk < r ^e iq < r ^/ 2 A a/3 (fc, q), (5) 

k.q 

Q) ^ ^ ^ / /3q:;77 / ^ )(' 0 q/ 2 +fc / ,7 < 0q/2 —fc',7') ? 

fc',77' 

( 6 ) 

with the nrulti-orbital interaction matrix V a p- nl '(k,k') 
and V’a(r) = exp(7fe • r). 

The multi-band BdG equations are then expressed as 

Jdr 2 H(r 1 ,r 2 ) < p(r 2 ) = E < p(r 1 ). (7) 

With the use of the eigenvectors </>(jt), the mean-field 
BdG Hamiltonian (HJ is diagonalized. 


A. Model and assumptions 


Let us discuss the decoupled multi-band quasiclassical 
theory, which is appropriate in the conventional two-band 
models for MgB 2 [20|. We assume that all TV x TV matrices 
in the BdG Hamiltonian are diagonal in the “band” basis, 
expressed as 


-*Vi) 


(17) 

J a,8 

A (ri,r 2 ) 

~ A a (ri,r 2 )5 a p. 

Ctfi 

(18) 


This assumption is equivalent to that off-diagonal inter¬ 
band elements are neglected. In this case, the normal 
state Hamiltonian in momentum space is expressed as 


ll N (fc) 


/Ai(fc) 0 0 

0 0 

\ 0 0 A N (k) 


(19) 









B. Appropriate region of the decoupled 
Eilenberger equations 
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with the eigenvalues A i(k). In real space, the Fourier 
transformation on each band makes the band-diagonal 
Hamiltonian. The mean-field BdG Hamiltonian in 
Eq. H]) is rewritten as 


H = ^ E f dr 1 dr 2 'Jf a (ri)^H a (r 1 ,r 2 )'if a (r 2 ), (20) 

a ^ 


with 


H a (ri,r 2 ) = 

( H^(ri,-iVi)S(ri-r 2 ) A„(ri,r 2 ) 

V A*(r 2 ,r!) ~H^*(r 1 ,-iV 1 )S(r 1 - r 2 ) 

( 21 ) 

Here, the column vector ’®' Q ,(r) is composed of fermionic 
annihilation %p a and creation operators ip' a on the a 
band at the position r (a = l,--- ,N), \I ’ a {r) = 
(ip a (r),ipl l (r)) T . There is no inter-band effect in the 
Hamiltonian, since the BdG Hamiltonian H a is deter¬ 
mined on each band. The gap equations in Eq. © are 
rewritten as 

^ ^ Vg / y (fc? ^ ) i'$q/2+k',‘y'4*q/2— ^', 7 )• (^ 2 ) 

fc, 7 

In this approximation, the multi-band effects are in¬ 
cluded as the inter-band pairing interactions only in the 
pairing V^ 7 . Thus, the quasiclassical decoupled Eilen¬ 
berger equations are easily obtained as 

iv F , a V R g R (k F ,z) + [za z - A R {k F )a z ,g R {k F ,z)]_ = 0, 

(23) 

A“(fc F ) = -^EE f ^VMk F ,k' F )f R (k' F ,iu; n ), 

n f3 J 1 F| 

(24) 


with quasiclassical Green’s function on the a band ex¬ 
pressed as 



Here, we neglect the self energies and the vector poten¬ 
tials for simplicity, z denotes a complex energy, er 2 is the 
Pauli matrix in 2 x 2 Nambu-Gor’kov space, [A, B\- = 
AB — BA is used, and G^(fc, z) is the Green’s function 
of the a band. It should be noted that the self-energy 
must be diagonalized in this basis and A in this section is 
2x2 matrix in Nambu space. The above equations were 
used in the simple multi-band superconductors such as 
MgB 2 . For MgB 2 , there are two bands (so-called er- and 
7r-bands) and the gap equations connect information on 
each band. 


Now, we discuss the appropriate region of the decou¬ 
pled Eilenberger equations in the previous section. We 
introduce the Hamiltonians FT orbltal (fc) and H hand (k) in 
normal states with the orbital basis and the band basis in 
momentum space, respectively. Generally, the Hamilto¬ 
nian H orhltal (k) has the off-diagonal elements, since the 
orbital basis is a basis which is diagonal in real space. 
With the use of the momentum dependent unitary ma¬ 
trix U(k), one can obtain the diagonal Hamiltonian ex¬ 
pressed as 



/ Ai(fc) 0 

° \ 


(k)H OThital (k)U(k) = 

0 

° 

, (27) 


^00 

Ajv(fe) / 


= H hand {k). 


(28) 


The “band” basis is determined by the unitary transfor¬ 
mation of the orbital basis in momentum space. If the 
unitary matrix does not depend on momentum, one can 
choose the basis which simultaneously diagonalizes the 
Hamiltonian in both real and momentum spaces. Gener¬ 
ally, the decoupled Eilenberger equations are derived by 
neglecting the off-diagonal elements. 

We show three examples that the decoupled Eilen¬ 
berger theory is not appropriate as follows. The first 
example is an impurity problem in the multi-band su¬ 
perconductors. In the previous study [71!), they assumed 
that the impurity-induced self energy was described as 
a momentum-independent band-diagonal matrix, which 
lead to the decoupled Eilenberger equations. We point 
out that this assumption induces non-local impurities in 
real space, as follows. In the band basis with this as¬ 
sumption, the Gor’kov equations become 

G h k and (z) = G°^ and {z) + G^ band (z)E band (z)G^ and (z). 

(29) 

Here, A band denotes the matrix defined by the ba¬ 
sis which diagonalizes the normal-state Hamiltonian 
.fforbital. To describe impurities in real space, one has 
to use the orbital basis in real space. With use of the 
unitary transformation from the band basis to the or¬ 
bital basis, the impurity-induced self-energy in the or¬ 
bital basis E orbltal should have a momentum dependence 
expressed as 

E orbital (fc, z) = U{k)t hand {z)U\k), (30) 

with 

If the self-energy in the band basis is obtained by the 
T-matrix approximation for randomly distributed point 
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impurities given as 

S band (z) = n imp V + n imp ^ VG° p oMtal (z)V , (32) 

V 

the self-energy in the orbital basis becomes 
£ orbital (fc, z) = n imp U orbital (fc, k) 

+ n imp F orbital (fc,p)G° orbital (^)^° rbital (p, k), (33) 

V 

with the effective “non-local” impurity potential defined 
as 

V olhiu \k,p) = U(k)VU\p). (34) 

Therefore, the decoupled Eilenberger equations does not 
describe the local impurity potentials. 

The second example is the proximity-induced 
impurity-robust p- wave effective order parameter on 
the surface of a topological insulator, as discussed 
later in Sec. IVII El With the use of the band basis, an 
effective chiral p-wave order parameter can be derived 
by the previous quasiclassical framework. This previous 
framework, however, can not describe the robustness 
against non-magnetic impurities, which was proposed by 
directly solving the BdG equationsjUj. The impurity 
robust p-wave superconductor is naturally introduced in 
our framework. 

The third example is the appearance condition of 
the Andreev bound states at a surface. In a single 
band model, the Andreev bound states occur when the 
sign of the gap function changes through the scattering 
process[z3j. In a multi band model, an ambiguity of the 
” sign” of the order parameter makes the above statement 
unclear. This ambiguity can not be resolved by the pre¬ 
vious quasiclassical framework. In our multi-band qua¬ 
siclassical Eilenberger approach, we can overcome this 
difficulty by deriving the most appropriate effective or¬ 
der parameter, which obeys the statement of the Andreev 
bound states as discussed later in Sec. EH 

We can use the decoupled Eilenberger equations if 
we assume the momentum-independent unitary matrix 
(U(k) = U ). In this assumption, however, we can not 
treat the “orbital” character. Therefore, we propose the 
general multi-band Eilenberger theory. 


IV. QUASICLASSICAL TREATMENT I: 

WAVE-FUNCTION APPROACH 

In this section, we derive the quasiclassical equations 
on the basis of the BdG equations. The quasiclassical 
theory is founded on an assumption that the coherence 
length £ is much longer than the Fermi wave length l/k F 
(i.e. t;k F "C 1)@- This assumption is valid, if the or¬ 
der parameter amplitude is much smaller than the Fermi 
energy, and this condition is fully fulfilled in BCS weak- 
coupling superconductivity. In this theory, the wave 
function is expressed by a product of the fast oscillating 


one characterized by the Fermi momentum p F and the 
slowly varying one by the coherence length. We proposed 
the quasiclassical theory for the multi-orbital topological 
superconductor fiij . The generalization of this theory is 
proposed in this section. 


A. Assumptions 


We assume that the eigen vector (f>(r) in Eq. J7]) is 
expressed by a product of the fast oscillating one charac¬ 
terized by the Fermi momentum and the slowly varying 
one by the coherence length expressed as 

<? H r ) = eIfcFr< ^fcp ( r )’ ( 35 ) 


where 


M , 

^kA r ) = J2 ( 

1=1 '■ 


uf(k F )ff*(r)\ 
v i N ( fc F)fff P (T-) ) ' 


(36) 


Here, / ; fcp (r) and gf p (r) correspond to slowly varying 
components, uf(k F ), v^(k F ) are the fast oscillating 
functions adopted as normal-state uniform eigenvectors 
satisfying the eigen-equations, 


H m (k) 


<(fc) 

«r(*o 



where 


H m (k ) 


H m (k) 0 \ 

0 H m (-k)* ) ' 


(37) 


(38) 


The Fermi surfaces in normal states are expressed by the 
set of the zero-energy eigenvalues of H N1 (k). The M- 
eigenvalues e/(fc) cross the Fermi level (i. e. ez(fcp) = 0). 
We assume that the eigenvalues near the Fermi level are 
same (ei(fc), ■ • • , £jf(fc) = £(fc)) expressed as 

* N1(fe) ( $$) =« fc) ($($)■ (39) 

This assumption is appropriate when M denotes the in¬ 
ternal degrees of freedom at the Fermi level as shown in 
Fig.ETa). We note that there is an exception as shown 
in Fig. [2Kb). For example, this exception occurs when 
the Fermi level is located at the center of the Dirac cone. 
However, the assumption is usually appropriate in many 
realistic materials. We should note that there is the re¬ 
lation between uf(k) and u ; N (fc) at the Fermi energy 
expressed as 


v?(k F )=u?*(-k F ). (40) 


B. Andreev-type equations 
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FIG. 2. The schematic figures of the electron bands in the 
multi-orbital system characterized by M = 2 at the Fermi 
wave momentum k f. 


By substituting Ea. (l35l) into the BdG equations (El), we obtain the Andreev-type quasiclassical equations. Eventu¬ 
ally, we have the 2 M x 2 M quasiclassical BdG equations represented as (in detail, see Appendix lAl). 


/ -iv F • V + V 0 (r,k F ) A eS (r,k F ) \ / / fcp (r) \ _ / / fcp (r) \ 

l AI ff (r,fe F ) ^ F -V-F 0 *(r,-fe F ); 7 \ 9 kF (r) ) ’ 


(41) 


with introducing M x M matrices Vq and A e fj defined by 


V 0 (r,k F ) = U^H m (r)U^, 
A eS (r,k F ) = U^A(r,k F )U%, 

(42) 

(43) 

with the N x M matrix given by 


U% = (u»(k F ),...,u* M (k F )). 

(44) 

Here, the M-component vectors f kp (r) and 
note f kF (r ) T = (/^(r), • • • ,/^ p (r)) and 
(ffi P ( r )) • • • : 9m ( r ))> respectively. The N x 
has the relation expressed as 

g kp (r) de- 
g fcF (r)T = 

M matrix 

U k u k — IjWxM' 

(45) 

Note that U k f U k f< ^ Inxn if M is not equal to N. With 
the use of the N x M matrix : the eigenvector 4>{r) 

in the BdG equations El is approximated as 

Mr) ~ T 1 

' A, U"V»(r) ) 

(46) 


The resultant quasiclassical BdG equations m are 
equivalent to the linearized BdG (Andreev) equations if 
we consider the single band system (N = 1) [9| . Thus, 
we successfully reduce the matrix dimension from the 
number of the bands N to the number of the degenerated 
Fermi levels M in this quasiclassical treatment. 


C. Boundary condition at a specular surface 

Let us discuss the boundary condition at a specular 
surface. For simplicity, we consider that the material is 


filled in the region z > 0. By assuming the translational 
symmetry along the x and y axes which conserves the 
momentum fe F || = (k Fx ,k Fy ), the boundary condition is 
given by 

<fr(k n ,z = 0) = 0. (47) 

First, we find the solutions which satisfy the above 
boundary condition in normal states at the Fermi energy, 
expressed as 

w N (fc F ||, z = 0) = 0. (48) 

By solving the eigenvalue equations with the normal- 

state N x N Hamiltonian H N1 (k): 

H m (k n ,ki)u^ kFhH) =0, (49) 

the boundary condition becomes 
K M 

c i( fc F||)<( fcp|| ,fci) = 0. (50) 

i l 


Note that k * l z is a complex number and satisfies Im k l z > 
0. Here, K denotes the number of k l z with the same 
conserved momentum (fc F ||). For example, in the case of 
K = 2 and M = 1, the boundary condition is satisfied 
only when the vector fe ^ is parallel to the vector 

u i,(k n ,k z ) expressed as 


u 


N 

2,(fe F || > k z) 


= e ^ 12 


u 


N 

l,(fe F ||,fc*)‘ 


(51) 


Here $12 is the overall phase difference between these 
two vectors. Thus, we obtain 

c 2 = -e- i4>12 ci. (52) 
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In a superconducting state, we use the quasiclassically- 
approximated wavefunctions. In the quasiclassical treat¬ 
ment, the eigenvector 0(fe F ||,r) in the BdG equations is 
approximated as 


<fr{k F \\,r) 


K 

Ye ik > 


l 


U% n , kl) f k k k F\\,z) \ 

U^ ki) 9 k *(k n ,z))’ 


with the boundary conditions 
K M 

EE c '( fe FII )/^( fc Fil ^ = o) = 0 ’ 

i l 
K M 

J2J2 cl ^~ k n)9^ z {k n ,z = 0 ) = 0 . 

i l 

Here, we define the N x M matrix given by 

tjM _ ( n ... ..n \ 

U {kw\\,K) ~ i M l,(fcF|| ,**)> 


(53) 


(54) 

(55) 


(56) 


In order to find the coefficients c((fc F ||), we have to solve 
the boundary condition (15U1) in normal states. For ex¬ 
ample, in the case of K = 2, the boundary condition in 
Eq. (15U1) is expressed as 


TT1V1 M , TTM M __ n 

u (k F[[ ,kl) c (k n ,kl) + u (k n ,kl) c (k n ,kl) ~ 


(57) 


with c (kl h k%) = ( c K fc F||)," ■ ,cf (fc F ||)). We show the 
boundary conditions of two superconducting systems 
with K = 2, which can be easily obtained, as exam¬ 
ples. In the case of N = M, the solution of the boundary 
condition (1571) is 


the above boundary condition can be only used when 
the normal-state eigenvectors with different momenta k\ 
and k z are parallel to each other shown in Eq. ED, 
since we have to find the correct boundary condition for 
/ fc ^(fe F ||, z = 0) which satisfies Eq. (1531) . 

The characteristic momentum k\ is obtained by solving 
a normal-state eigenvalue equation (1491) with the bound¬ 
ary condition (1501) at the Fermi energy. This momentum 
k l z is usually a real number since wavefunctions with a 
Fermi wave number can usually satisfy the boundary con¬ 
dition. In this case, we solve the Andreev equations (TUI) 
with the real-number momentum (k F ^,kl). However, 
the momentum k l z is not the Fermi wave number, when 
there are bound states in normal states. The normal- 
state wavefunction localized at a boundary is described 
by a complex momentum. We will show the system which 
needs the complex wave number to describe the bound 
states, as discussed in Sec. IVIII 

We discuss the existence condition of solutions in 
Eq. (15U1) . With the use of the matrix representation, the 
equation (1501) is rewritten as 

UC = 0, (64) 

where the N x MK matrix U and MA'-dimension vector 
C are defined as 


U = (^(fc PHl fcj). • • ■ - ))’ 

( C M \ 


C = 


rM 


\ c (fe F „,fcf) J 


(65) 

( 66 ) 


M _/r M t tjM M 

C ( fc F|| ,kl) ^(fcpn ,fc|) C (fcpj| ,fcj) ’ 


(58) 


since the relation = Inxn is satisfied. 

Note that Eq. (1551) is not the solution in Eq. (1571) if N ^ 
M. By substituting Eq. (1551) into Eqs. (1531) and (1551) . we 
obtain 


f k Hz = 0) = -vg*;;$f k Hz = 0), (59) 

9 k Hz = 0) = -^ ( ( X7,uf *g k hz = 0), (60) 

where the “transfer matrix” V’// Sp|| ’ fc 1 V is determined as 

(fepn 


T>(feF||,*£) _ frM\ f T M 

(fcp|| (kF\\,k%) (fepih^i)' 


(61) 


Next, we consider the case of M = 1 and K = 2. The 
boundary condition (1551) becomes 


The above equations are linear homogeneous equations 
with MK unknowns. These equations can have the so¬ 
lutions when MK > TV. 

Finally, we note that the general boundary condition 
for the conventional quasiclassical equations has been dis¬ 
cussed by several groups [H,@. The incoming and out¬ 
going wavefunctions are connected by the scattering ma¬ 
trix S expressed as 

0(&out) = *Sfc out fe in (/>(fci n ), (67) 

where fc in ( out ) is the wave number of the incoming (out¬ 
going) quasiparticles. Here, S is the scattering matrix 
defined at the Fermi energy in normal states[64|. The 
above relation can be used to develop the general bound¬ 
ary condition in the multi-band superconductors. The 
development of the general boundary condition for the 
multi-band quasiclassical equations is a future issue. 


f k *(k F \\,z = 0) = -e l * 12 f k *{k F §,z = 0), (62) 

9 k Hk n ,z = 0) = -e^ 2 g k *(k Fh z = 0), (63) 


V. QUASICLASSICAL TREATMENT II: 
GREEN’S FUNCTION APPROACH 


which is equivalent to the boundary condition in the sin- In this section, we derive the equation of motion of 
gle band case when $12 = 0. We should note again that the quasiclassical Green’s function so-called Eilenberger 




equations in a multi-band system. The size of the matri¬ 
ces of the Green’s function is reduced by the low-energy 
projection. 


We note that there is another Gor’kov equation called the 
“right-hand Gor’kov” equation. In terms of the Wigner 
representation, the right-hand equation is expressed as 


A. Wigner representation 


The Wigner representation is usually introduced in 
terms of the derivation of the quantum Boltzmann equa¬ 
tions. The transport-like equations of motions of the qua- 
siclassical Green’s functions in the single band system are 
systematically derived with the use of the Wigner rep¬ 
resentations. We introduce the Wigner representation 
defined by 


A{R, k) 


dfe~ lk ^A{r u r 2 ) 


ri—R-\-r/2,r2—R—v/2 


( 68 ) 


Here, R = (ri + r 2 )/2 and r = r± — r 2 are the center-of- 
mass coordinate and the relative coordinate, respectively. 

The Gor’kov equations in the Wigner representation 
are expressed by 

(ioj n - H no (R) - H m {k ) - A (R, k ) - t(R, fc, iuj n )) 

* G(R, k, iuj n ) = 1. (69) 

Here, we introduce the ^-product (Moyal product) deter¬ 
mined by 


A(R , k) * B(R, k) = exp 

x A(R,k)B(R',k') 


- (V*/ -Vfl - Vfc ■ Vjj') 


k'=k.R'=R 


• (70) 


G(i£, fc, iu n )-k 

(iu n - H m (R) - H m (k) - A(R, k ) - E (R, fc, iu n )) = I. 

(71) 


The local density of states with the Wigner representa¬ 
tion is expressed as 


N(r,E) = —Im 

7r 


lim 


•n-t+o 




Tr G(r , k 1 iu> n —> E + irj) 


(72) 


Let us derive the quasiclassical equations from Eq. m 
as follows. In superconductors, the characteristic length 
of the center-of-mass coordinates is the coherence length 
£, which is much longer than that of the relative coor¬ 
dinates characterized by 1/fcF- Assuming that the char¬ 
acteristic coherence length is long, the Moyal product in 
the first order of Vr is given as 


A{R, k) * B(R, k) ~ A{R, k)B{R , k) 

+ l - (V r A(R, k) ■ V k B{R, k) - V k A{R, k) • V r B{R, k)) . 

(73) 

Then, the Gor’kov equations are expressed as (see, Ap¬ 
pendix m 


(■ ico n ~ H m {R) - H m (fc) - A {R, k) - t{R, k, iu n )) G(R, k, ico n ) 

- l -V R [H m (R) + A (R, k) + S(R, k, ito n )] ■ V k G(R , k, iu n ) + ^ V k H m (k) ■ V r G(R , k, iu n ) = 1. (74) 


The above equations are simultaneous differential equa¬ 
tions with k and R. In a single-band system, the above 
equations becomes the differential equations with respect 
to R with a parameter fcp (i.e. the quasiclassical Eilen- 
berger equations), since we can eliminate the fc-mixing 
term VfcG(J?, k, iu> n ), with the use of the contour inte¬ 
gration with respect to the energy and the relation 
Vfc oc d/d^k- Thus, in order to derive the multi-band 
quasiclassical Eilenberger equations, the fc-mixing term 
V k G{R,k,iuj n ) has to be eliminated. However, we can 
not simply integrate the above equations with respect to 
As, since At is the energy obtained by diagonalizing the 
Hamiltonian H N1 (k) and the Gor’kov equations are de¬ 
termined in the orbital basis. One has to characterize the 


equation by the Fermi wave momentum in the low-energy 
region. 


B. Projection to the effective low-energy model 

We introduce the projection matrix to develop the 
multi-band quasiclassical theory. The projection matrix 
eliminates the degree of freedom about the high energy 
region. We define the 2 N x 2 N projection matrix as 

P k = Gf G fc Mt . 


(75) 
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Here, the 2N x 2M matrix U k is defined by 


G(k,ioj n ) — P k G(k,iuj n ) becomes 


t4 M = 


u£f o \ 

0 U M * ) ’ 


(76) 


where the N x M matrix is defined by Eq. (Iddl) . The 
projection operator satisfies the relation 

P k P k = P k , (77) 


since the matrix U k ! always satisfies the relation 


[6G(k,iOJn)\ 


N-M 


E 


Ucn(k)U^(k) 
iu n — e 7 (fc) 


(82) 


with 1 < a, f3 < N. If the eigenvalues are located far 
from the Fermi energy (|w n | -C e 7 (fcp)), 5G(k F ,ioJ n ) be¬ 
comes negligible small so that we can obtain the relation 
expressed as 


v?'v? 


= i 


2Mx2M- 


(78) 


P k G(k , iu> n ) ~ G(fc, iu n ), (83) 


In order to show the physical meaning of the projec¬ 
tion matrix, we operate P k on the homogeneous iV-band 
Green’s function in normal states determined by 

G{k,iu n ) = (iuj n i-H m (k))- 1 . (79) 


The component of G(k,iui n ) is expressed as 


U a Jk)m(k) 

G a0 (k,iuj n ) = ^2 —— - t,.\ ’ (1 <«,/?< IV) 


rt iw n -e 7 (fc) 


(80) 


Here, the TV x TV unitary matrix U(k) diagonalizes 
T7 N 1 (fc) and e 7 (/c) is the y-th eigenvalue of TL N 1 (fc). Op¬ 
erating Pfc on G(k,iu) n ), we obtain 


which is appropriate in the low-energy region. 

We introduce the 2Mx2M reduced matrix A expressed 
as 


A(k) = U?'A{k)U*, (84) 

where A is the 2N x 2TV matrix used in the Gor’kov 
equation (f^Ul) . With the use of Eqs. (1751) and (1551) . A can 
be expressed by 


A{k) ~ U k A{k)U k \ (85) 

in the low-energy region (fc ~ k F ). 

C. Multi-band quasiclassical Green’s function 


[P k G(k, iuj n )\ a p 


E 


u m 

^ hex 7 


frM* 

u k(3j 

-ffc 


(81) 


with 1 < a, (3 < TV. Here, we use the assumption in 
(1551) . The above equation means that the projection 
operator Pk eliminates the information of large eigen¬ 
values from G(k,i<jJ n ). The difference 5G(k,iuj n ) = 

I 


Let us construct the 2 M x 2 M quasiclassical multi¬ 
band Eilenberger equations in the projected space. By 
multiplying the the both sides in Eq. (El by the matrices 
U k and U k , subtracting the right-hand Gor’kov equa¬ 
tion, and integrating over (in detail, see Appendix [Cl) . 
we obtain 2 M x 2 M quasiclassical multi-band Eilenberger 
equation expressed as 


iv F (k F ) • V R g R (k F ,z ) + \za z - V 0R {k F )cr z - A R {k F )a z - T, R (k F , z)a z , g R {k F , z)\ _ = 0, ( 86 ) 


where we introduce the 2 M x 2 M Green’s function, non¬ 
local potentials, order parameters, and self-energies de¬ 
termined by 


G R (k,z) = U^G(R,k,z)U f 

(87) 

V 0R (k) = U^H m (R)U™, 

( 88 ) 

A R (k) = U^A(R,k)U^, 

(89) 

E R (k,z) = U^t(R,k,z)U^. 

(90) 


Here, cr z denotes the Pauli matrix in the Nambu-Gor’kov 
space and we define the 2 M x 2 M quasiclassical Green’s 


function expressed as 

g R {k F ,z) = j) d£ k a z G(R,k,z), (91) 

where j> means taking the contributions from poles close 
to the Fermi surface. The matrix structure of the above 
equation is equivalent to that of the single-band Eilen¬ 
berger equation. When the eigenvalue is not degenerated 
(M = 1), the 2x2 equation (l 86 l) can be regarded as 
that in a spin-singlet single band superconductor [46M48| . 
Therefore, we call the 2 M x 2 M matrix representation 
the “single-band description”. The band index a is useful 
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to compare with the present multi-band theory by replac¬ 
ing g^(k F ,z) as g R (k F ,z). The multi-band Eilenberger 
equations (1551) are similar to the decoupled Eilenberger 
equations Eq. (1551) . We should note that our multi-band 
theory includes the orbital characters, since all matrices 
are determined in the projected space. For example, the 
self-energy with the T-matrix approximation depends on 
the momentum as discussed in Eq. (|82)>. 

We should note that the normalization condition is 
equivalent to that in the single-band system expressed 
as (in detail, see Appendix iDl) 

gg=-n 2 l. (92) 


D. Relations in the projected space 

Let us discuss the relations satisfied in the pro¬ 
jected space. M x M order parameter is defined in 
Eq. (1431) . If the original order parameter have the re¬ 
lation A (R, fe) 1 'A(i?, fe) = A 0 (J?, k)l NxN , the projected 
order parameter have the similar relation: 

A e fi(R, ky A e g(R, k) ~ A 0 (R, k)l M xM, (93) 

because of U^Up' A{R, k) ~ A (R, k) near the Fermi 
energy. Here, A e g(R,k) is determined in Eq. (1451) . This 
indicates that the unitarity of the order parameter is con¬ 
served in the projected space. 


E. Physical quantities and Gap equations 


Let us express physical quantities with the use of the 
multi-band quasiclassical Green’s function. By substitut¬ 
ing Eq. (1551) into Eq. m, the local density of states is 
expressed as 


N(r, E) = —-Im 

7T 


— 1 

= —Im 

7r 


lim V Tr G 11 (r,k, iw n —> E + ig) , 

77—»+0 zJ 

k 

(94) 

)( lim<Tr g n (r, k, iu n -> E + ig)) kw , 

(95) 


where G 11 and g 11 denote (1, l)-element 
in the particle-hole space, the bracket de¬ 
notes the Fermi-surface average (• ■ ■ )f e = 

/ • • ■ d5 , F(fcF)|’UF(fcF)| -1 / J d5 *F(^f)|^f(^f)| —1 ■ Other 
physical quantities can be expressed by the multi-band 
quasiclassical Green’s function in the same manner. 

Finally, we complete the multi-band quasiclassical the¬ 
ory by giving the self-consistent equations for the order 
parameters. The gap equations in the Wigner represen¬ 


tation is given as 

N 

A a p(R,k) = -T y y Vpa-j-y' (fc, k’)F iy (R, k ', iu n ) 

fc',7 7 ' n 

(96) 

By using Eq. (l89|) . the M x M order parameter matrix 
is expressed by 


AeffiduH^F) = 

M 

- ( 0? ) 
n Z 3 Z 4 

where g 12 denotes (1,2)-element in the particle-hole 
space and (k F , k’ F ) is the effective interaction ex¬ 
pressed as 


y^ 2 (fc F ,4) = Vp a] yff(kF, & f ) 


ot.fi ll' 


frM* frM* 
X U k F ah U 


tM* 


frM tti\ 

— k F pi 2 u k' F -fl 3 u — fcp7'i 4 ■ 


(98) 


We can simplify the above gap equations if the pairing 
interaction has a separable form expressed as[65j 


V Parrr (k,k') = V a p{k)V ryf {k > ). (99) 

Here, we use the relation V^ a;77 '(fc, fc') = 
L-y 7 ;a/'3 (k , fe)E3. The gap equations in this case 
are given as 


M 

Ar = (^F )9l3UR (100) 

n l 3 l A 

with 

A e s R (k F ) = A R V(k F ), (101) 

V{k F ) = U^V(k F )U%, (102) 

where we assume that Vyy(k')* = E 77 /(fc'). 


F. Perturbative approach in the quasiclassical 
theory: Zeeman and spin-orbit couplings 

In this section, we discuss the method to treat the Zee- 
man and spin-orbit couplings. In the previous studies [3^. 
[35l . f66l468i |, they used the 4 x 4 matrix quasiclassical Eilen¬ 
berger equations in spin and Nambu spaces expressed as 


iv F ■ V R g R {k Fl z) 

+ \za z - Act z - Hi(k F )a z ,g R (k F ,z)\_ = 0. (103) 

Here, cq denotes the Pauli matrix in the Nambu space 
and Hi(k F ) includes the spin-orbit and/or Zeeman cou¬ 
pling terms. These quasiclassical equations can describe 
a vortex state in a Fulde-Ferrell-Larkin-Ovchinnikov 
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superconductor [681. In terms of our theory, the above 
equations are obtained by assuming that the number of 
degenerated Fermi surfaces is two (M = 2). In general, 
however, the Zeeman and spin-orbit interactions split the 
degenerated bands (i. e. M = 2 —*■ M = 1 ). Thus, it is 
found that the previous studies assume that the Zeeman 
and spin-orbit interactions are weak. With the use of 
the perturbative approach in the multi-band quasiclassi- 
cal theory, we can derive the above equations as follows. 
Let us divide H N 1 (fc) in Eq. (l39ll into the two terms as 

H m (k) = H$\k) + Hi(k). (104) 

Thus, in order to construct the projection operator P &, 
we can use the eigenvectors obtained by 

4 N1 (fcK = & 4 - (105) 

This approach is appropriate for the case that the inter¬ 
band pairing between the different Fermi surfaces is im¬ 
portant. The perturbative Zeeman held enables us to 
treat the Pauli paramagnetic depairing in the quasiclas- 
sical framework. 


G. Riccati-type equations 


straight line in the direction of the Fermi velocity vp: 

vf— = — 2u n a — aA^a + A. (HO) 

Os 

fij, _ _ 

vp — = 2oj n b + bAb — A^. (Ill) 

os 

In a bulk system with A^A oc 1, the solutions of the 
Riccati equations are 

d(w n ) = - A ==, (112) 

Wn + \J^n + ^Tr AAt 

At 

b(u n ) =- . ==. (113) 

w « + V w n+ ^Tr AAt 

According to the previous paper[l2j, putting a(s a ) = 0 
and b(sb ) = 0 as initial values, one can obtain physical 
solutions by integrating the Riccati equations. Here, s a 
and Sb are initial spatial points. 


H. Boundary condition for the Riccati parameters 
at a specular surface 


It is known that it is difficult to numerically solve the 
quasiclassical Eilenberger equationsfUJ, since the equa¬ 
tions have a divergent solution as a particular solution. 
A careful computational treatment is required for inte¬ 
grating the Eilenberger equations with the use of the so- 
called explosion method [69(. To avoid this difficulty, the 
Riccati-type equations, which are obtained by a special 
parametrization form of the quasiclassical Green’s func¬ 
tion, are used[37i-[39l. [4ll-[43j. In addition, to solve the 
Riccati equation stably, we have proposed the efficient 
numerical method for obtaining unique solutions in the 
single-band Eilenberger framework[l2j. We show that it 
is easy to expand this method into the multi-band sys¬ 
tems. For simplicity, we neglect the self-energy E = 0. 
We use a special parametrization form of the quasiclas¬ 
sical Green’s function to solve Eq. ([ 86 ]) . The solution g 
of Eq. (l 86 l) can be written as, 


g = —ittN 


(1 — ab) 2ia \ 
— 2 ib — (1 — ba ) ) ’ 


N = 


(l+d6) 1 0 

0 (l + 6d) -1 


(106) 

(107) 


where M x M matrices a and b are the solutions of the 
following matrix-type Riccati differential equations: 


To solve Riccati-type differential equations, one has to 
consider the boundary condition for Riccati parameters 
a and b. In this section, we consider a specular surface. 
In the case of M = 1 or M = 2, we can show the trans¬ 
formation of the linearized BdG equations to the matrix 
Riccati equations. Thus, the boundary condition for the 
Riccati parameters can be determined explicitly. We note 
that, in many materials, the number of the degenerated 
Fermi levels M is not larger than M = 2. 

If the Fermi level with a certain momentum fcp is not 
degenerated (M = 1), the Riccati equations are derived 
by the relation: 




(114) 


With the use of the boundary condition for wavefunctions 
in Eqs. (l62l) and (l63l) . the boundary condition at a surface 
z = 0 with K = 2 is given by 


~a{kl z ) = e~ 2i ^ 12 a(kp z ). 


(115) 


If the Fermi level is doubly degenerated (M = 2), we find 
the relation expressed as 


a(r, &f) = iU(kp, r)V (A:f, r) _ 


(116) 


up • Va = —2 uj n d — aA^a + A. (108) 

vp ■ V 6 = 2uj n b + bAb — A^. (109) 

Since the above equations contain V only through up • V, 
these can be reduced to a one-dimensional problem on a 


where the 2x2 matrices U(kp,r ) and V{kp 1 r) are de¬ 
termined by 


U(kp,r) = ( / fcp (r) -g fcp (r)* ) , (117) 

V(kp,r) = (g kF (r) f~ kF (r)*). (118) 
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We obtain the boundary condition at a surface z = 0 
with K = 2 and N = M = 2 expressed as 


a \ K Fz) ~ V {k Fa ,,k Fv ,k\, z ) a \ K Fz) 


(^Fi i^Fy ’^Fz ) 


(119) 


Note that T4 


(^Fi ?&F y 5^Fz ) 


(^Fa; y i^p 2 ) 


is determined in Eq. (ED). 


I. Arbitrary transformation about normal-state 
eigenvectors: Appearance condition of the Andreev 
bound states 


We discuss an arbitrariness of the N x M matrix U kF . 
With the use of this arbitrariness, one can discuss the 
appearance condition of the Andreev bound states in 
multi-band superconductors. As shown in Eq. (Iddl) . the 
N x M matrix consists of the zero-energy degener¬ 
ated eigen vectors about the normal-state Hamiltonian 
H N1 (k). Because of the degeneracy, it should be noted 
that the N x M matrix U kF has additional degrees of 
freedom expressed as 

< = ^ 4 , ( 120 ) 


with a M x M arbitrary unitary matrix A kp . Although 
this matrix does not change the 2N x 2 M projection 
matrix P kp , the representation of the effective order pa¬ 
rameters A e ff(r,fcF) determined in Eq. (l43l) depends on 
the matrix A kF , expressed as 

A eff (r, k F ) = A{ F U^A(r, k F )U%A*_ kp . (121) 

It should be noted that this arbitrary transformation 
does not change any physical quantities, since the ma¬ 
trix Ak F changes the boundary condition. With the use 
of the unitary matrix A kp , we can simplify the boundary 
condition as follows. 

In the case of K = 2 and N = M, the boundary con¬ 
ditions (15U1) and (16U1) become 




where 


T ?,( fe F||,^) _ i T>(feFH,*2) if 

v (fcp||,fci) - W fc F|i AD '(fepn,fcj) A {k n M z )' 

By using the matrix Ak F satisfying the relation 
^(*=fh Af) = ^(fcF||AI)(^(fc P || ll ,fci) ) 


( 122 ) 


n ' z> g k *(z = 0), (123) 


(124) 


(fc Fii- fe 2)M (125) 


we obtain the simplified boundary condition expressed as 


f k l( z = 0) = -f k *(z = 0), (126) 

g *Z(z = 0) = -g*{z = 0). (127) 


In addition, in the case of K = M = N = 2, the bound¬ 
ary condition for the Riccati amplitudes (111911 becomes 

a(k Fz ) = a(k^ z ), (128) 

The above boundary condition is equivalent to that for 
spin-triplet superconductivity in the past quasiclassical 
treatment. In the case of an effective one-band system 
(M = 1) with K = 2 , the lxl unitary matrix A kp is 
rewritten as A kp = e l ^ k v . Thus, we can erase the overall 
phase 4 >i 2 in Eq. (l52l) . The boundary condition becomes 

f k *( z = 0) = -f k *{z = 0), (129) 

g k *(z = 0) = -g k *(z = 0). (130) 

We obtain the boundary condition for the Riccati ampli¬ 
tude expressed as 

a(k Fz ) = a(fc£j, (131) 


which is completely equivalent to the boundary condition 
in the single-band quasiclassical Eilenberger theory. 

Now, we can discuss the appearance condition of the 
Andreev bound state at a surface in multi-band super¬ 
conductors. In a single band model, the Andreev bound 
states occur when the sign of the gap function changes 
through the scattering process. In the case of K = 2 
and M = 1, we can easily discuss the appearance condi¬ 
tion of the Andreev bound states. Note that the bound 
states appear if the condition (l5ll) is satisfied in this case. 
To use the above boundary conditions (11311) . the order 
parameter matrix after the scattering process should be 


A e ff(r, k F ) 





(132) 


with 



- fr M h~r M 
= u ki u ky 


(133) 


The quasiclassical Green’s function at a surface diverges 
when the relation 


1 + a(ito n — > e + irj, kp)b(ioj n —)■ e + irj, fc|) = 0 (134) 

is satisfied. 0 With the use of the bulk solutions in 
Eqs. (11121) and (11131) . the appearance condition of the 
zero-energy Andreev bound states becomes 

|A off (fci)||Aeff(fe|)| + A off (fci )A off (fe|)* = 0. (135) 

Thus, the sign of the order parameter is important for the 
appearance condition of the Andreev bound state even in 
multi-band superconductors. 


VI. MULTI-BAND EFFECTS 

We discuss the physical meanings of the multi-band 
quasiclassical theory described by Eq. (15S1) . In our theory, 
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FIG. 3. The schematic figure of the electron bands in the 
three-orbital superconductors. 


the “multi-band” effect is characterized by two factors. 
The first factor is how many solutions are in Eq. (1591) 
(i.e. the information of the eigenvalues). The second 
one is how the orbitals are mixed in Eq. (1591) (i.e. the 
information of the eigenvectors). 


A. Eigenvalues 

We discuss the eigenvalues obtained in Eq. (l39ll . In the 
quasiclassical theory, the number of the solutions at the 
Fermi surface M characterizes the multi-band effect. For 
example, in the three-orbital spin-singlet superconduc¬ 
tors, the superconducting order parameter is described 
by a 3 x 3 matrix (N = 3) in Eq. ©. Let us consider 
that the only one band crosses the Fermi energy at the 
Fermi surface with the Fermi wave number k F (M = 1) 
as shown in Fig. [3] In this case, the other bands must 
have the much higher or lower energies. Thus, the super¬ 
conducting order parameters on or between these bands 
(e. g. An (fcp) or Ai 2 (fc F )) can not affects the physi¬ 
cal quantities, since these order parameters (~ meV) are 
much smaller than the energy scales of electron bands (~ 
eV). Therefore, the lxl order parameter matrix on the 
band crossing the Fermi energy (e. g. A 22 (fc F ) in Fig. [3]) 
is only effective for the superconductivity. In terms of the 
above point, many multi-band superconductors such as 
MgB 2 or iron-pnictides can be described by the single¬ 
band superconducting gap because of M = 1 in these 
materials. 


B. Eigenvectors 

We discuss the eigenvectors in Eq. (1551) . which de¬ 
scribes the ratio of the hybridization of the orbital char¬ 
acters at the Fermi wave number fcp. The multi-band 
effects are described by these eigenvectors. For example, 
we consider the impurity self-energy with the Born ap¬ 
proximation. In our framework, 2 M x 2 M self-energy 


becomes 

T,(k F ,iuj n )a z = ni m pV (fcp, k F ) 

+ n imp (V (fc F , k F )g(k Fl iu n )V (k F , k F )) k > F , (136) 

where we determine 

V(k F ,k' F ) = U^V 0 a z U^. (137) 

In the case of an effective single band system (M = 
1), V(k F ,k F ) becomes the 2x2 matrix and its (1,1)- 
component is expressed as 

F n (fcp, fcp) = u\(k F )V 0 ui(k F ) (138) 

with 

i7 N 1 (fcp)ui(fcp) = 0. (139) 

In the case of Vo oc 1, the strength of the multi-band ef¬ 
fects can be determined by the orthogonality between the 
eigenvectors at the different Fermi momenta. We note 
the case of the T-matrix approximation for randomly dis¬ 
tributed impurities. The 2 N x 2 N matrix self-energy is 
written as 

E (iuj n )a z = n imp f(iu> n ), (140) 

where 

T{ico n ) = Vo + Y, Gw (iu n )f(iu n ). (141) 

k' 

In our framework, the self-energy with 2 M x 2 M matrix 
form is obtained as 

E(/cp,ica n ) —n.i m pXfc; F fc F (zcc n ), (142) 

where 

Tk F k’ F (iwn) = V (fcpi fcp) + {V (fcp, k F )g(k F , iu>n)Tk^.,k' F )k^ ■ 

(143) 

Thus, the eigenvectors of the normal state Hamiltonian 
are important to describe the impurity effects. 

C. Non-local anisotropic potentials in the 
projected space 

As shown in the previous sections, the multi-band 
Eilenberger equations (l86l) have the orbital characters 
through the matrix U M (k F ). It should be noted that the 
matrix U M (k F ) makes the potential H N0 (R) non-local 
and anisotropic. The non-locality originates from the 
non-unitary transformation which erases the information. 
Non-local potentials have been used as pseudo-potentials 
in the first-principles calculations. In the first-principle 
calculations, the pseudo-potential method which treats 
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valence electrons only is commonly used to erase the de¬ 
grees of inner-shell electrons. In our theory, the non¬ 
local potentials are used to erase the fast oscillations 
characterized by kp. The multi-band effects are under¬ 
stood by the non-locality and anisotropy in the projected 
space. The effective potential Von.(kp) in Eq. (l86l) is non¬ 
local and anisotropic, since the potential depends on the 
center-of-mass coordinate R and the relative coordinate 
kp. The non-locality and anisotropy can be easily un¬ 
derstood through the example of the self-energy with the 
Born equation expressed as Eq. (I136D . The above self¬ 
energy can be regarded as that made from the non-local 
potential V(r,r'). We should note that the non-locality 
and anisotropy are strong in the iron-based superconduc¬ 
tors, since d-orbitals are strongly entangled at the Fermi 
level. 


D. Differences between the present and previous 
quasiclassical multi-band frameworks 

We discuss the differences between the present and pre¬ 
vious quasiclassical multi-band treatments. Fundamen¬ 
tally, note that our framework is an extension of a previ¬ 
ous single-band Eilenberger framework. Thus, we make 
the derivation similar to the previous single-band one. 

Our main point is the low-energy projection which sys¬ 
tematically reduces a JV-band system to a M-band sys¬ 
tem. There are two kinds of the reductions in our paper. 
The first one is same as that in the previous paper [47f, 
which uses Fermi velocities and Green’s functions at the 
Fermi level in normal states. The second one is the re¬ 
duction of the band-degree of freedom, which was not 
mentioned in Ref. et al. Our equations can treat the off- 
diagonal elements of a Green’s function between bands. 
A previous theory treated inter-band effects only through 
gap-equations not the Eilenberger equations. Note that 
the results in the previous papers using the quasiclassical 
theory were obtained only in the problem which the pre¬ 
vious framework was available. Our framework extends 
the applicable region of the quasiclassical theory. 

We claim that there was no multi-band Eilenberger 
equation which can treat inter-band effects correctly in 
inhomogeneous systems. In the previous framework, for 
example, it is hard to study the vortex bound states with 
impurities in complicated multi-band superconductors, 
such as iron-based superconductors. One usually consid¬ 
ers the five-orbital model for the iron-based superconduc¬ 
tors whose Fermi surfaces are constructed by only three- 
bands. In this case, the decoupled Eilenberger equations 
have three band indices. When one introduces a self 
energy ( e.g. the impurity-induced self energy), the off- 
diagonal elements of the self energy can not be described 
by the decoupled Eilenberger equations. On the other 
hand, such a self-energy matrix should be defined by 
the Dyson equation with Feynman-diagram techniques 
in five-orbital model. In addition, band-coupled Eilen¬ 
berger equations except for our framework can be con¬ 


structed only in the case that the Fermi velocities are the 
same in the different bands, since the decoupled Eilen¬ 
berger equations with a band index are characterized by 
a Fermi velocity on the band. 

Finally, we show the example which makes clear differ¬ 
ence between our and previous frameworks, qualitatively. 
The corrections in our framework describe the multi¬ 
orbital effects. Multi-band effects in the present frame¬ 
work are necessary to correctly calculate a reduction of 
the critical temperature caused by impurities. While the 
self energy with T-matrix approximation induced by the 
randomly-distributed impurities does not depend on the 
momentum in the previous decoupled Eilenberger theory, 
the self-energy in our Eilenberger equation depends on 
momentum. If the momentum dependence is neglected, 
the quasiclassical theory can not correctly calculate the 
reduction of the critical temperature. We show the 
qualitative difference between the present and previous 
frameworks in the system on the topological insulator in 
Sec. lVII El The theoretical calculation by directly solving 
the BdG equations suggested that the proximity-induced 
superconductivity on the surface of the topological insu¬ 
lator is robust against non-magnetic impurities The 
present multi-band framework can correctly describe the 
robustness. The previous quasiclassical framework, how¬ 
ever, can not reproduce this robustness. 

VII. MULTI-BAND QUASICLASSICAL 

APPROXIMATIONS IN VARIOUS KINDS OF 
SYSTEMS: EXAMPLES 

In this section, we apply the multi-band quasiclassical 
theory in the various kinds of systems as examples. 

A. Noncentrosymmetric Superconductors: CePtsSi 

We show that our multi-band quasiclassical theory 
makes the past debates clear. The noncentrosymmetric 
superconductor CePtaSi has the Rashba-type spin-orbit 
interaction due to the lack of the inversion symmetry [TuL 
l49j . The mixed spin-singlet-triplet model has been used 
to study this material. By assuming that the spatial 
variations of the s-wave pairing component of the pair 
potential are the same as those of the p-wave pairing 
component, the gap function is expressed as 

A(fc) = [4'er 0 + d k ■ er] icr y . (144) 

Here, cq is the Pauli matrix in spin space. The normal- 
state Hamiltonian with the Rashba-type spin-orbit cou¬ 
pling is written as 

H(k) = (X k - n)a 0 +g k ■ cr, (145) 

where the spin-orbit interaction satisfies the relation 
= — £/&. Here, A*, is the dispersion without the spin- 
orbit interaction. We determine gk as [Tol l49j 

g £ = g{— sin ^ sin 0, cos ^ sin 0,0). (146) 
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We assume that the d-vector is parallel to g k expressed 
as 


d k = Ad{— sin sin 0, cos 0 sin 0,0). (147) 

In the previous papers [ 0 ], from the original Eilcnberger 
equation for noncentrosymmetric superconductivityI50L 
l53j , they have obtained two equations corresponding to 
these split Fermi surface I and II in the case of the s+p- 
wave pairing state, M 

*ui,ii ■ Vgpn + [iui n f 3 - Aijpgiji] = 0. (148) 


where A UI = [(n + if 2 )A hU - (fi - if 2 )A I * II ]/2, A I>n = 
ip± Ad sin 9 are the order parameters on the Fermi surface 
I and II, (fi, 72 , 73 ) are the Pauli matrices in the particle- 
hole space, and the commutator [a, 6 ] = ab — ha. There 
are many successes with the use of the above decoupled 
equations. We should note that there are some debates 
about the appropriate region of the above approach [52} . 
In the real material such as CePtsSi, there is the strong 
spin-orbit coupling (~ eV). It has been not clear whether 
this approach is the weak-spin-orbit coupling approach, 
since the two same-size spherical Fermi surfaces are as¬ 
sumed. The difference of the size of the each Fermi sur¬ 
face depends on the strength of the spin-orbit coupling. 
On the other hand, the size of the Fermi surfaces is nat¬ 
urally considered in our multi-band theory. 

Let us apply the multi-band quasiclassical theory to 
the noncentrosymmetric superconductors in order to de¬ 
rive the decoupled Eilenberger equations Eq. (I148D di¬ 
rectly from the Hamiltonian (11451) . The eigenvalues of 
the 2x2 matrix in Eq. (11451) is given by 

t±(k) = \ k — p±\g k \. (149) 


Although there are two Fermi surfaces, the eigenvalue is 
not degenerated so that we obtain M = 1. The eigenvec¬ 
tors associated with e±(fc) are expressed as 


u t = 


1 

V2 

1 

72 



(150) 

(151) 


The lxl effective gap function is given as 

A± = ±ie ±i<f, {^± A d sin(9). (152) 

The above effective gap function is not equivalent to Apn 
in the previous paper. We should note that a representa¬ 
tion of the effective gap function has an arbitrary degree 
of freedom expressed as 


A c ff(0, 9)' = 9) f A cff ((/), 9)A((f> + 7r,9 + 7r)*, (153) 


as discussed in Sec. IVII Thus, we can use a 1 x 1 arbitrary 
unitary matrix in order to change a representation of the 
effective gap. With the use of the lxl unitary matrix 
A± (</>, 9) defined as 


±2 


1 — 


the effective gap function can be rewritten as 


A' ± = $ ± A d sin 9, (155) 

which is equivalent to that in the previous papers. 0,0 
Next, we assume the degenerated Fermi surface (M = 
N = 2), which is appropriate when \gk\ -C 1. With the 
use of the unitary matrix A(</), 9) defined as 


A{^9) 



(156) 


the effective gap function can be rewritten as 


A eff(tM)' 


f + A d sin 9 

0 


0 

T - A 


d sine 


(157) 


In terms of the multi-band quasiclassical theory, we clar¬ 
ify that the decoupled equations (11481) are valid with the 
arbitrary strength spin-orbit coupling. 

Finally, we consider a specular reflection at a surface 
perpendicular to x — y plane. We consider that the quasi¬ 
particles before and after a scattering have momentum 
fci = ( (j>i,9), fc 2 = {<j>2,9), respectively. By assuming 
the degenerated Fermi surface (M = N = K = 2), the 
unitary matrix with A(</>, 9) in Eq. (I156D becomes 




M 


1 

72 


-ie-^l 2 -e "^/ 2 \ 
e ^/2 ie i<t ,/2 J - 


(158) 


whose effective order parameter is given in Eq. (11571) . 
The transfer matrix VJ *“ 2 = is expressed as 


t>fc 2 _ f cos A</> — sin A(j> \ 
fc i l sinA</> cos A 4> J ’ 


(159) 


with A cj) = {<f>\ — (f> 2 )/2. This transfer matrix suggests 
that both intra- and inter-band scatterings occur at a 
specular surface. The surface bound states and spin cur¬ 
rents discussed in the previous paper 0 ) can be explained 
in terms of this band-active surface. 


B. Three-orbital model: S^RuCL 

Let us apply our theory to a multi-band supercon¬ 
ductor. In this section, we consider Sr 2 Ru 04 as the 
three-band superconductor. The many tight-binding 
models for Sr 2 Ru 04 have been proposed by several 
authors0, [53457} . According to Ref. 3f| the effec¬ 
tive tight-binding Hamiltonian is expressed by the three- 
orbital model characterized by d yz ~, d xz -, and d xy - or¬ 
bitals (N = 3) expressed as 

/ t v k z -n ef+iX -A \ 

H(k) = ef - iX e x k * — g iX , 

V -A -*A e^-p) 


A±((j>,9) = ±e 


(154) 


(160) 
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FIG. 4. (Color online) Fermi surfaces for the three-band su¬ 
perconductor Si' 2 Ru 04 . 


FIG. 5. (Color online) The in-plane anisotropy of the effective 
potential 9')\ at the most inner Fermi surface in the 
three-band superconductor Sr 2 Ru 04 . 


where 

e v ^ = —2 t 2 cos (k x ) — 2fi cos (k y ), 
e%. z = —2ti cos (k x ) — 2 1 2 co s(k y ), 
e k V = -2t 3 (cos{h x ) + cos (k y )) - 4t 4 cos (k x ) cos (k y ) 

— 2 * 5 (cos( 2 fc x ) + cos (2 k y )), 

= —4*6sin(fc x ) sin(fcj / ), (161) 

with A = 0.032, t[ = 0.145, t 2 = 0.016, t 3 = 0.081, f 4 = 
0.039, *5 = 0.005, *6 = 0, and p, = 0.122. Here, we adopt 
the material parameters in Ref. [36, which can successfully 
describe the three Fermi surfaces for Sr 2 Ru 04 as shown 
in FigQ] We call each band as band I, band II, and band 
III in ascending order of the eigenvalues. 

We consider the non-magnetic impurities to discuss the 
non-local anisotropic effective potential. We introduce 
the in-plane anisotropy of the effective potential defined 
as 


these systems, the layer-dependent spin-orbit coupling 
induces the exotic superconducting states. In the N- 
layer spin-singlet s-wave superconductor, the multi-band 
Eilenberger equations with 4iV x 41V matrix quasiclassical 
Green’s function are written as 


iv f • V<?+ [zt z - A- K(k F ),g\_ = 0, 


with 


K(k F ) — ^*_L-^inter 4“ z^ *0 

(H so (k F ) 0 

V 0 H* s 0 (-M J ’ 


+ 


A = 


po®Inxn 0 

0 ~o'o®Inxn 

0 A 

-At 0 


(163) 


(164) 

(165) 

(166) 


V(M') = «L p(fl )«M«')> ( 162 ) 

Here, k F (d) denotes the position of the most inner Fermi 
surface (i.e. band III) in momentum space (k F (0) = 
k F (8)(cos 9, sin 9)). As shown in Fig. 0 the right-angled 
scatterings suppress in the most inner Fermi surface. 
This suppression originates from the fact that the eigen¬ 
vector associated with k F (8 = 0 ) mainly consists of d xz 
orbital and the eigenvector associated with k F (8' = 0 ) 
mainly consists of d yz orbital. 

C. Heavy fermion CeCoIns/YbCoIns superlattice: 

The perturbative approach 

In this section, we consider the system with both the 
spin-orbit coupling and the Zeeman interaction. The 
locally noncentrosymmetric systems are realized in the 
heavy fermion CeCoIn 5 /YbCoIn 5 superlattice[ 6 (J. In 


.ffinter = CTO ® T±, (167) 

Hz = —o z ®Inxn, (168) 

Hso(k F ) = g(k F ) ■ a <g> S d . (169) 

Here, a z is the 2x2 Pauli matrix, Inxn is the unit 
matrix, T± is the hopping matrix between layers, and 
Sd = diag(aj, • • • , ay) denotes the layer-dependent spin- 
orbit interaction. In the above equations, the hopping, 
the Zeeman, and the spin-orbit coupling terms are re¬ 
garded as the perturbations with 21V-degenerated Fermi 
surfaces. With the use of this perturbation theory, one 
can treat the inhomogeneous system with vortices [6 lj. 

D. Topological superconductors with the strong 
spin-orbit coupling: Cu^B^Ses 

We discuss the boundary condition in the three- 
dimensional topological superconductor with the strong 
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spin-orbit coupling in this section. CuzB^Sea is the 
one of the candidates of the topological superconductors 
where the topologically protected Majorana bound states 
form at the boundary. We have proposed the quasi- 
classical framework on topological superconductors with 
strong spin-orbit coupling [il]. In the previous paper jjj), 
we have obtained the linearized BdG equations called An¬ 
dreev equations by decomposing the slow varying com¬ 
ponent from the total quasi-particle wave function. Ap¬ 
plying this quasiclassical treatment, the original massive 
Dirac BdG Hamiltonian derived from the tight-binding 
model represented by 8 x 8 matrix is reduced to 4 x 4 
one. The resultant Andreev equations become equiva¬ 
lent to those of spin-singlet or triplet superconductors 
without the spin-orbit coupling. In this section, we show 
that the same result is obtained by the Green’s function 
techniques. 

The normal-states effective Hamiltonian for Cu x Bi 2 Se 3 
is expressed as 

Here cq denotes the Pauli matrix in the spin space. In 
the quasiclassical theory, H{k) is regarded as H N1 (k) 
in Eq. ®. The eigenvalues of the 4x4 matrix are 
degenerated expressed as 


ei(fe) = ±E 0 {k) - /x, (171) 


consider M(k) = Mo(k\\) + M\k 2 . By using Eq. (l49l) . we 
find that the wave numbers become 

( ^ ±)2 = 2^[ 1 + 2MoMl±a 1 ]’ (176) 


with £ fc|| = y/(l + 2M 0 M 1 ) 2 + 4M 0 (g 2 - - |fc|| | 2 ). 

When the condition g > yj— |fc||| 2 is satisfied, there 
are two real wave numbers and two imaginary wave num¬ 
ber, expressed as 

k\_ = k _ (177) 

(178) 

k\+ = ir)+ , (179) 

k 2 z+ = - ir )+, (180) 


with k- 



1 + 2AfoMi-kJ /(2 M 2 ) 


and rj + = 


yj [l + 2MqMi + /(2 M 2 ). By assuming that the 

material is filled in the region £ > 0, the coefficient of the 
wavefunction with k 2 , is zero. Thus, we obtain K = 3. 
The boundary condition (l50l) is expressed as 


( M (fc||A-) M (fc||,-fc-) U (k ll ,ir 1+ ) ) c 

+ ( u \k h k_) u (fc lh -fc_) u tk v z V+ ) ) C 2 = 0. (181) 


where Eo(k) = \JM{k) 2 + \k\ 2 . In the case of fj > 0, 
the eigenvectors u k and u k (M = 2), and eigenvalue £& 
in Eq. (15^1) are given as 


Zk = E 0 (k) - g 


(172) 


u l k = c 


Xi 

kcr 

E 0 (k)+M(k ) 



(173) 


where xT = (1,0), xi = (0,1), c = ^/( E o + M(k))/2E 0 . 
We consider the 4x4 odd-parity fully-gapped gap func¬ 
tion (so-called A lu state) defined as 


A s A ° (4 7 ) ■ t 174 > 

By substituting this gap function into Eq. (l43l) . we obtain 
the 2x2 effective gap function written as 


When we consider fc|| = 0, the coefficients are given as 


c^c 2 


(* 

( 1 1 M+M(fe_) IT)- 

1 ~r /i+M(+ + ) 

)\ 

1 

/ fj.+M(k _) ir)- 
/i+M(+ + ) fc_ 


2 


V 

-1 

/ 


(182) 


Thus, the boundary condition for the quasiclassical wave 
function is obtained as 


f~ k ~ (z = 0) = cif k ~ (z = 0), 
g ~ k - ( z = 0 ) = c \g k - (z = 0 ), 

where 

_ k_(g + M(ir/ + )) - ig+(g + M(k-)) 
Cl fc_(/i + M(iri + )) + ir] + (g + M(fc_)) ’ 


(183) 

(184) 


(185) 


A ° ff(fc ) = l^(k) k ' <T (* <Ty )' (I 75 ) 

This gap function is completely equivalent to that in the 
previous paper [44j in terms of the Dirac BdG Hamilto¬ 
nian. 

Let us consider the boundary condition with a specular 
surface at z = 0. We adopt the boundary condition that 
all components of the wave-function becomes zero at z = 
0, which is different from that discussed in Ref. [3l|. We 


In the non-relativistic limit ( \k- \ <C M(fc_)), the bound¬ 
ary condition becomes 

f~ k -{ z = 0)=-f k -(z = 0), (186) 

9 k ~ (z = 0) = —g k ~ {z = 0), (187) 

which is equivalent to that in a single band quasiclassical 
framework. This result is consistent with the fact that 
this superconductor becomes a p-wave superconductor in 
the non-relativistic limit (32|. 
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E. Robust p-wave superconductivity on a surface of 
topological insulator 


In this section, we discuss the impurity effect in the 
s-wave gap superconductivity on a surface of topological 
insulator. Let us show that the proximity-induced su¬ 
perconductivity on a surface of topological insulator is 
robust against nonmagnetic impurities [ 63 . We consider 
that an s-wave superconductor is deposited on the sur¬ 
face of the topological insulator [70j. The effective two- 
dimensional Hamiltonian on the surface is described as 


with 


H(k) 


ho(k) A itjy \ 

A*(-ia y ) -h%(-k) ) 


(188) 


ho(k) = vk ■ <7 — /ioo- (189) 


The eigenvalues of the 2x2 normal state Hamiltonian 
h 0 (k) are given by 

e±(k) = — n ± v\k\. (190) 

The eigenvectors associated with e±(k) are expressed as 


< = I®) 

= yf () ’ ( 192 ( 


where k = ( k x , k y ) = k( cos <j>, sin <j>). In the case of /r > 0, 
the lxl effective gap is given as 

A eff (fc) = Ae*. (193) 


Thus, the effective p- wave superconductivity appears on 
a surface of the topological insulator. 

Let us consider a nonmagnetic impurity effect in 
this proximity-induced superconductor. The non- 
perturbative quasiclassical Green function in a homoge¬ 
neous system is given as 


g(k F ,ioj n ) = 


iuj n Ae^ 
y/w* + |A| 2 V -Ae-^ -iu. 


. (194) 


The effective potential V(k F ,k' F ) in Eq. (113711 is ex¬ 
pressed as 


V(k F ,k’ F ) = V 


p i<50/2 


cos 


(f) 


0 

e -i54>/2 cog 


mr 


(195) 


with 5(f> = (j> — <j>'. Here, V is the amplitude of the po¬ 
tentials. The second order of the impurity self-energy in 
Eq. (11361) becomes 


r 2 

£(fc F , A; n ) (2) CT z = n imp V 2 g(k F ,iu] n ) 

Jo 


dcj)' cos 2 


27T 


(196) 


Since this self-energy satisfies 

[E(fepj iwn)^ 2 *^, g(k F , ioj n )\- = 0, the quasiclassi¬ 

cal Eilenberger equations with the self-energy are 
completely equivalent to those without the self-energy. 
Therefore, this proximity-induced superconductivity on 
the surface of a topological insulator is robust against 
nonmagnetic impurities. 

Finally, we point out that the previous decoupled 
quasiclassical framework can not reproduce the robust¬ 
ness against nonmagnetic impurities proposed in Ref. 1(52] . 
With the use of the band basis, the effective gap is given 
in Eq. (11931) at the Fermi energy. Since this effective gap 
means p- wave superconductivity, the superconductivity 
should be fragile against nonmagnetic impurities in the 
previous decoupled quasiclassical framework. 


F. Surface quasiclassical theory: the partial 
quasiclassical approximation for topological 
insulators 

Let us discuss the “partial” quasiclassical approxima¬ 
tion with considering the topological insulators. The su¬ 
perconductivity in surface states on topological insulators 
has been attracted much attention because of the stage of 
the Majorana Fermion and a quantum computing. With 
the use of the proximity effects from the superconductor 
on the topological insulator, the two-dimensional mass¬ 
less Dirac quasiparticles due to the surface bound states 
on the topological insulator form the superconducting 
Cooper pairs. Therefore, it is important to construct the 
two-dimensional effective Eilenberger equations originat¬ 
ing from the normal-state surface bound states. 

Let us consider the surface at z = 0. By introducing 
the coordinate r = ( x,y,z ) = (r±,z), we can define the 
partial Wigner representation expressed as 

^ZiZ2 (-^_L ; k^_ ) = 

J dr j_ e~ ik± ' f± - A ZlZ2 (r ± + . (197) 


Here, R± = (ru + r 2 j_)/2 and r± = ru — r 2 j_ are the 
center-of-mass coordinate and the relative coordinate on 
the two-dimensional plane parallel to the surface, respec¬ 
tively. The projection operator is determined as 

P k±ZlZ2 =U^( Zl )U^(z 2 ), (198) 

where 



with the matrix U k ^(z) = (ul ± (z), ■ ■ ■ ,u^f ± (z)). Here, 
the vector u l k (z) is the eigenvector expressed as 

I? N1 (fcx,zK» = Z^uljz). (200) 
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We note that the Hamiltonian H N1 (k,z) includes infor¬ 
mation about a presence of a surface. The projected 
effective gap function is given by 

A (R±,k ± ) = j dz^U^iz^A^R^k^U^^), 

( 201 ) 

since the projection includes the z-integration written as 

Azi2 2 = J dz 3 Pk ZlZ 3 A Z3Z2 . (202) 


Let us consider the three dimensional topological super¬ 
conductor as an example. The eigenvector in Eq. (11701) 
with the boundary condition u l k (z = 0) = 0 is expressed 
as mi 


(*) 


e 2M i sinh(ATz) 

2 VA 


( e-^ \ 



(203) 


where £ k± = ^jk% + k% - /x, fcj_ = (. k x ,k y ) = 

,Jk% + ky (cos 4>, sin </>), M{k) = Mo(k±) + Mik z , 

K = (Vl + 4M 0 Mi)/(2Mi), and A 
f 0 dzexp(z/M 1 )\ sinh(A"z)| 2 . Here, we assume p > 0 
and obtain the above solution with the use of the 
perturbation with respect to k±. By substituting the 
eigenvector in Eq. (12031) and the odd-parity fully-gapped 
gap function in Eq. (11741) into Eq. (12011) . we obtain 


A(J?j.,fej.) = 0. (204) 


Therefore, the proximity-induced odd-parity fully- 
gapped gap function does not open the spectral gap as 
shown in Ref. [33l 

With the use of the above method, we directly show 
that the proximity-induced s-wave superconductivity on 
the topological insulator can be regarded as a chiral p- 
wave superconductivity, as shown in Sec. IVII El By sub¬ 
stituting the eigenvector in Eq. (12031) and the even-parity 
fully-gapped spin-singlet intra-orbital gap function ex¬ 
pressed as 

A s A “ (7 £)' < 205 > 

into Eq. (12011) . we obtain 

A(R ± ,k ± ) = Aj_e # , (206) 

which is equivalent to the chiral p -wave superconductiv¬ 
ity in Eq. (11931) . 


G. Others 

Finally, we discuss the advantage of the multi-band 
quasiclassical theory. The computational cost drastically 


decreases with the use of our theory in multi-band sys¬ 
tems. Thus, we can treat the inhomogeneous systems 
such as those with vortices and surfaces easily, in order 
to discuss the magnetic field dependence of the multi¬ 
band superconductors. Since we do not use any assump¬ 
tions about the electronic structures in normal states, the 
superconducting system with the arbitrary tight-binding 
Hamiltonian derived by the first-principle calculation can 
be mapped onto the effective low-energy system. In the 
theoretical point of view, one might develop the general 
theory for impurity effects in multi-band superconduc¬ 
tors, since the multi band effects are explicitly included 
as the non-local and anisotropic potentials. It should 
be noted that one can understand what is neglected in 
the quasiclassical theory in multi-band superconductors. 
One might know the difference between the single-band 
and multi-band superconductors through the study with 
our multi-band Eilenberger theory. 


VIII. SUMMARY 

In summary, we proposed the unified quasiclassical 
multi-band Eilenberger equations in order to map the 
multi-band systems onto the effective systems in the re¬ 
duced space. We derived both the Andreev and Eilen¬ 
berger equations with an arbitrary boundary condition. 
We showed that the resultant multi-band Eilenberger 
equations are similar to the single-band ones, except for 
some corrections to describe multi-band effects. The or¬ 
bital characters on the Fermi surfaces in normal states 
are included in our theory. Our theory could describe the 
past studies with the use of the quasiclassical Eilenberger 
theory. Since we do not use any assumptions about the 
electronic structures in normal states, the superconduct¬ 
ing system with the arbitrary tight-binding Hamiltonian 
derived by the first-principle calculation can be mapped 
onto the effective low-energy system. The potentials, 
the order-parameters, and self-energies in multi-band sys¬ 
tems were mapped onto the non-local ones in the reduced 
space as shown in Eqs. (|87|) - ([90l) . We showed that the 
self-energy with the T-matrix approximation of the non¬ 
magnetic impurities becomes non-local and anisotropic. 
We pointed out that this non-locality is similar to the 
pseudo potential in the first-principles calculations. The 
multi-band effects can be understood by the non-locality 
and the anisotropy in the mapped systems. 
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Appendix A: Derivation of the Andreev equations 


We derive the multi-band Andreev-type equations as 
follows. We substitute Ea. (l35l) into the BdG equations 

I 

The diagonal blocks in the BdG equations become 




ikpri nkp 


fr(n)u?(k F ) 


ik^r i 


tf wo (ri)/f p (ri)-i 


dH N1 (k) 


dk 


V./f p (n) 


k—kj? 


«"(**)• (Al) 


The off-diagonal blocks are converted as 


dr 2 A(r 1 ,r 2 ) 


e ik ^g^(r 2 )vf(k F ) 


e ik ^A{r u k F )9^{ri)vf{k F ), 


(A2) 


where A(r l5 k F ) is the order parameter matrix in the Wigner representation. 

I- 


With the use of the relation 

dm 


r?M t dH N1 (k) ~ M 
Uk dk k 


fc=fe F dk 


k—kp 


L MxM 


= fpl 


MxM) 


(A3) 


we obtain the multi-band Andreev-type equations (Idll) . 


Appendix B: Expansion of the ^-products 

We expand the ^-products for A (R, k) and E(.R, k) 
written as 

[A (R, k) + E (R, k)] * G(R , fc, iu n ) 


In the quasiclassical theory of superconductivity, the *- 
products for H N0 (R) and i7 N1 (fe) are expanded in the 
1 st order written as 

H m {R) * G(R, k, iuj n ) ~ H m (R)G{R, k, ioj n ) 

+ *-V R H m (R) ■ V k G(R,k,iuj n ), 

<B2) 

H m (k)*G(R,k,iu n ) ~ H N1 (k)G(R, k, iuj n ) 

- l -V k H m (k)- V R G(R,k,iu n ). 

(B3) 


Appendix C: Derivation of the multi-band 
quasiclassical Eilenberger equations 


r _ „ , x We derive the multi-band quasiclassical Eilenberger 

[&(R,k) + X(R,k)]G(R,k,iu n ). (Bl) equations . 


By multiplying the the both sides in Eq. m by the matrices U and Ujf , we obtain 


(iu n - V 0 {R, k) - £ k a z - A (JR, k) - E (R, k , iu n j) G{R, k , iu n ) 

- ^V R V 0 (R,k) ■ If*' [V k G(R,k,iuj n )] Ujf + l -U^ [V k H N1 (k)] U fe M • V R G(R,k,ico n ) = I. 


(Cl) 


The right-hand Gor’kov equation in the projected space are written as 

G(R, k, iu n ) ( iuj n - V 0 {R, k) - £ k d z - A (JR, k) - E(JR, k , iu n )) 


-Gf f [V k G(Rk,ioj n )\ U™ ■ V R V 0 (R,k) --V R G(Rk,ioj n ) ■ [V k H m (k)] Ujf = I. 


(C2) 
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By subtracting the right-hand Gor’kov equation, the equation (IC1I) become 

[icj n a z - H 0 (.R, k)a z - A(R, k)a z - S(J?, k)a z ,a z G(R, k; iu n )\ _ 


1 

2 L 


V R V 0 (R,k)a z ,a z U k 


m t 


[VfeG(J2, fc, ic< 


U, 


M 


+ + 2 L 


t/fc f [V fe i? N 1 (fc)] U^a z ,a z V R G(Rk-iu n ) 


= 0 , 


(C3) 


where [A, .B]± = AR ± BA. In the quasiclassical theory, the information on the 

Fermi surface is most important. By assuming that 
V 0 (R,k), A(R,k), E(i?, k), and U k are slowly varying 
functions around the Fermi momentum k F , we obtain 


[iUnVz - V 0 (R, k F )a z - A(R, k F )a z - E(FJ, k F )a zl a z G(R, k; iu n )] _ 


2 L 


V R V 0 (R,k F )a z ,a z U^ [V k G(R, k, iw„)] 


M 


J + 


- iv(k F ) ■ fc; iw„) = 0. 


(C4) 


The ^-integration erases the term with V k in Eq. (IC4I) , 
which is same procedure in the single-band model. Thus, 
the resultant 2 M x 2 M quasiclassical multi-band Eilen- 
berger equation is 


Appendix D: Normalization condition 

We consider the normalization condition for g. When 
g satisfies Eq. (1561) . g' = a I + (3g and dA are also the 
solutions of Eq. (l 86 l) as shown in Ref. [HJ. Thus, the 
solution of the equation has the form 


gg = al+/3g. (Dl) 


The coefficients a and ft are determined in a homoge¬ 
neous case. The 2 N x 2 N matrix Green’s function in a 
homogeneous system is written as [45] 


G(k,i<jj n ) 


UR) 


A + (fc,*w n ) B{k,iu n ) 

B*{k,iw n ) A_{k,iu n ) 


U\k), 

(D2) 


where 

A±(k,iuj n ) 

B(k,iu> n ) 


a/3 

a/3 


_ c i € a 

~ “^|A Q |2 + ( JW „)2_ e 2’ 

_ , A a 

a/3 -|A a |2 + ( JW „)2_ e 2. 


(D3) 

(D4) 


Here, the 2 N x 27V matrix U{k) is the unitary matrix 
which diagonalizes H N0 (k). Assuming that intraband 


pairing are dominant, we neglect the off-diagonal (inter¬ 
band) elements of the order-parameter matrix. 

The 2 M x 2 M Green’s function in the projected space 
G(k,iu} n ) is written as 


G(k,iu n ) 


A+ (k, iui n ) B M (k,ioj n )\ 
B M ^(k, ioj n ) AR(k,iUn)J' 


(D5) 


where 


A± (k : icj n ) 

B M (kGu n ) 


- a(3 


= s n 


iuj n ± £k 


= 6 , 


-|A a \ 2 + {iUn? - a 
A« 


a/3" 


a/3 — | Aq, | 2 + (iw ra ) 2 — 


2 ' 


(D 6 ) 

(D7) 


By substituting the above Green’s function into Eq. m, 
we obtain the multi-band quasiclassical Green’s function 
in a homogeneous system expressed as 


g(k F ,ioj n ) 

where 


a^(k F ,iu n ) b M {k F ,iu n ) \ , . 

-5 Mt (fc F ,iw n ) -a^(k F ,iu n ) ) ’ 


\a±{k F ,iu} n )\ a 

b M (k F ,iu n ) 


3 a/37T 


-6 t 


a/3 7T' 


%/^ + IAal 2 

A a 

V^n + |A a | 2 


(D9) 

(DIO) 


Therefore, the normalization condition becomes Ea. (l92l) 
even in an inhomogeneous system. 


[1] K. Kobayashi, M. Machida, Y. Ota, and F. Nori, Phys. [2] H. Kontani, and S. Onari, Phys. Rev. Lett. 104, 157001 
Rev. B 88, 224517 (2013) and references therein. 






























22 


( 2010 ). 

[3] K. Kuroki, S. Onari, R. Arita, H. usui, Y. Tanaka, H. 
Kontani, and H. Aoki, Phys. Rev. Lett. 101 , 087004 
(2008). 

[4] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 

( 2010 ). 

[5] Y. Ando, J. Phys. Soc. Jpn. 82 102001 (2013). 

[6] A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Lud¬ 
wig, Phys. Rev. B 78, 195125 (2008). 

[7] L. Fu and E. Berg, Phys. Rev. Lett. 105, 097001 (2010). 

[8] Y. Nagai, N. Hayashi, N. Nakai, H. Nakamura, and M. 
Okumura, and M. Machida, New J. Phys., 10, 103026 
(2008). 

[9] G. E. Volovik: Pisma Zh.Eksp.Teor.Fiz. 70 (1999) 601; 
JETP Lett. 70 (1999) 609. 

[10] Y. Nagai, Y. Ueno, Y. Kato and N. Hayashi: J. Phys. 
Soc. Jpn. 75 (2006) 104701. 

[11] G. Eilenberger: Z. Phys. 214 (1968) 195. 

[12] Y. Nagai, K. Tanaka, and N. Hayashi: Phys. Rev. B 86 
(2012) 094526. 

[13] P. Miranovic, M. Ichioka, and K. Machida: Phys. Rev. B 
70 (2004) 104510. 

[14] A. S. Melnikov, D. A. Ryzhov, and M. A. Silaev: Phys. 
Rev. B 78 (2008) 064513. 

[15] Y. Nagai and N. Hayashi: Phys. Rev. Lett. 101 (2008) 
097001. 

[16] Y. Nagai, Y. Kato and N. Hayashi: J. Phys. Soc. Jpn. 
75 (2006) 043706. 

[17] S. Graser, C. Iniotakis, T. Dahm, and N. Schopohl: Phys. 
Rev. Lett. 93 (2004) 247001. 

[18] C. Iniotakis, N. Hayashi, Y. Sawa, T. Yokoyama, U. 
May, Y. Tanaka, and M. Sigrist: Phys. Rev. B 76 (2007) 
012501. 

[19] N. Kopnin, Theory of Nonequilibrium Superconductivity. 
(Clarendon Press, Oxford, 2001). 

[20] M. Mugikura, H. Kusunose, and Y. Kuramoto, J. Phys. 
Soc. Jpn. 74, 1806 (2005). 

[21] K. Tanaka and M. Eschrig, Supercond. Sci. Technol. 22 
014001 (2009). 

[22] A Gumann, S Graser, T Dahm, and N Schopohl, Phys. 
Rev. B 73, 104506 (2006). 

[23] V. G. Kogan, C. Martin, and R. Prozorov, Phys. Rev. B 
80, 014507 (2009). 

[24] R. Prozorov, V.G. Kogan, Rep. Prog. Phys. 74, 124505 

( 2011 ). 

[25] M. Silaev, E. Babaev, Phys. Rev. B 84, 094515 (2011). 

[26] A. Moor, A. F. Volkov, and K. B. Efetov, Phys. Rev. B 
83, 134524 (2011). 

[27] A. F. Kemper, T. A. Maier, S. Graser, H.-P. Cheng, P. J. 
Hirschfeld and D. J. Scalapino, New J. Phys. 12, 073030 
( 2010 ). 

[28] M. Sato, Y. Takahashi, and S. Fuiimoto, Phys. Rev. B 

82, 134521 (2010) 

[29] Y. Nagai, Y. Ota, and M. Machida, J. Phys. Soc. Jpn. 

83, (2014) 094722. 

[30] M. Ishikado, Y. Nagai, K. Kodama, R. Kajimoto, M. 
Nakamura, Y. Inamura, S. Wakimoto, H. Nakamura, 
M. Machida, K. Suzuki, H. Usui, K. Kuroki, A. Iyo, 
H. Eisaki, M. Arai, and S. Shamoto, Phys. Rev. B 84, 
144517 (2011). 

[31] Y. Nagai, H. Nakamura, and M. Machida, J. Phys. Soc. 
Jpn. 83, 053705 (2014). 

[32] Y. Nagai, Y. Ota, and M. Machida, Phys. Rev. B 89, 
214506 (2014). 


[33] Y. Nagai, H. Nakamura, and M. Machida, JPS Conf. 
Proc. 3, 015013 (2014). 

[34] N. Hayashi, K. Wakabayashi, P. A. Frigeri, and M. 
Sigrist, Phys. Rev. B 73, 092508 (2006). 

[35] C. T. Rieck, K. Scharnberg, and N. Schopohl, .J. Low 
Temp. Phys. 84, 381 (1991). 

[36] V. B. Zabolotnyy, D. V. Evtushinsky, A. A. Ko- 
rdyuk, T. K. Kim, E. Carleschi, B. P. Doyle, R. Fit¬ 
tipaldi, M. Cuoco, A. Vecchione, and S. V. Borisenko, 
arXiv: 1212.3994 

[37] Y. Kato, J. Phys. Soc. Jpn. 69, 3378 (2000). 

[38] Y. Nagato, K. Nagai, and J. Hara, J. Low Temp. Phys. 
93, 33 (1993). 

[39] S. Higashitani and K. Nagai, J. Phys. Soc. Jpn. 64, 549 
(1995). 

[40] M. Sigrist and K. Ueda: Rev. Mod. Phys. 63 (1991) 239. 

[41] Y. Nagato, S. Higashitani, K. Yamada, and K. Nagai, J. 
Low Temp. Phys. 103, 1 (1996). 

[42] N. Schopohl and K. Maki, Phys. Rev. B 52, 490 (1995). 

[43] N. Schopohl, e-print arXiv:cond-mat/980464. 

[44] Y. Nagai, H. Nakamura, and M. Machida, J. Phys. Soc. 
Jpn. 83, (2014) 053705. 

[45] Y. Nagai and N. Hayashi, Phys. Rev. B 79, 224508 
(2009). 

[46] C. H. Choi and J. A. Sauls, Phys. Rev. B 48, 13684 
(1993). 

[47] J. W. Serene and D. Rainer, Phys. Rep. 101, 221 (1983). 

[48] Y. Nagai, Y. Ueno, Y. Kato and N. Hayashi, J. Phys. 
Soc. Jpn. 75, 104701 (2006). 

[49] N. Hayashi, Y. Kato, P. A. Frigeri, K. Wakabayashi, and 
M. Sigrist, Physica C 437-438, 96 (2006). 

[50] N. Hayashi, K. Wakabayashi, P. A. Frigeri, and M. 
Sigrist, Phys. Rev. B 73, 092508 (2006). 

[51] N. Hayashi, K. Wakabayashi, P. A. Frigeri, and M. 
Sigrist, Phys. Rev. B 73, 024504 (2006). 

[52] N. Hayashi, in private communication (2013). 

[53] N. Schopohl, J. Low Temp. Phys. 41, 409 (1980). 

[54] C. M. Puetter and H.-Y. Kee, Europhys. Lett. 98, 27010 

( 2012 ). 

[55] K. K. Ng and M. Sigrist, Europhys. Lett. 49, 473 (2000). 

[56] C. M. Puetter, J. G. Rau, and H.-Y. Kee, Phys. Rev. B 
81, 081105 (2010). 

[57] W.-C. Lee, D. P. Arovas, and C. Wu, Phys. Rev. B 81, 
184403 (2010). 

[58] B. de Wit and J. Smith, Field Theory in Particle Physics 
. (Elsevier, 1986), p.457 

[59] A.B. Vorontsov, I. Vekhter, M. Eschrig, Phys. Rev. Lett. 
101. 127003 (2008). 

[60] Y. Mizukami, H. Shishido, T. Shibauchi, M. Shimozawa, 
S. Yasumoto, D. Watanabe, M. Yamashita, H. Ikeda, T. 
Terashima, H. Kontani, Y. Matsuda, Nat. Phys. 7, 849 
( 2011 ). 

[61] Y Higashi, Y Nagai, T Yoshida and Y Yanase, J. Phys.: 
Conf. Ser. 568, 022018 (2014). 

[62] Y. Ito, Y. Yamaji, and M. Imada, J. Phys. Soc. .Jpn. 80, 
063704 (2011). 

[63] A. Shelankov and M. Ozana, Phys. Rev. B 61, 7077 

( 2000 ). 

[64] M. Eschrig, Phys. Rev. B 80, 134511 (2009). 

[65] P. B. Allen and B. Mitrovic, Solid State Phys. 37, 1 
(1983). 

[66] J.A.X. Alexander, T.P. Orlando, D. Rainer, and P.M. 
Tedrow, Phys. Rev. B 31, 5811 (1985). 

[67] A.B. Vorontsov and I. Vekhter, Phys. Rev. B 81, 094527 


23 


( 2010 ). 

[68] M. Ichioka, H. Adac.hi, T. Mizushima, and K. Machida, 
Phys. Rev. B 76, 014503 (2007). 

[691 E. Thuneberg, J. Kurkiiarvi, and D. Rainer, Phys. Rev. 
Lett. 48, 1853 (1982). 

[70] L. Fu and C.L. Kane, Phys. Rev. Lett. 100, 096407 


(2008). 

[71] P. Belova, M. Safonchik, K. B. Traito, E. Lahderanta, J. 
Phys. Conf. Ser., 303 012113 (2011). 

[72] S. Kashiwaya and Y. Tanaka, Rep. Prog. Phys. 63, 1641 

( 2000 ). 


